###############################################
#                                             #
# Value of Precision                          #
# Functions Necessary for Data Analysis       #
# Last Updated: 7/19/17                       #
#                                             #
# Josh Baker                                  #
# jbak@sas.upenn.edu                          #
#                                             #
###############################################


###########################################################################################
#                                   Instructions                                          #
###########################################################################################

## To load all functions in this document into your workspace, simply select the entire document
## and run all lines of code in R (or R-Studio).


###########################################################################################
#              Function for Identifying Which WEP Bin a Forecast Falls Within             #
###########################################################################################
id.wep <- function(val.unrounded){
  splits <- seq(from=0, to=1, 1/7)
  splits <- splits[-1]

  if(is.element(val.unrounded, splits)){
    val.unrounded <- val.unrounded + runif(1, min=-0.000001, max=0.000001)
  }
  
  wep <- ifelse(val.unrounded > splits[length(splits)], 7, min(which(splits > val.unrounded)))
  
  return(wep)
  
  #rm(val.unrounded,splits,wep)
}


###########################################################################################
#              Function for Identifying Which WEP Bin a Forecast Falls Within             #
###########################################################################################
id.nic <- function(val.unrounded){
  splits <- c(0.05,0.20,0.45,0.55,0.80,0.95)
  
  if(is.element(val.unrounded, splits)){
    val.unrounded <- val.unrounded + runif(1, min=-0.000001, max=0.000001)
  }
  
  nic <- ifelse(val.unrounded > 0.95, 7, min(which(splits > val.unrounded)))
  
  return(nic)
  
  #rm(val.unrounded,splits,nic)
}


###########################################################################################
#              Function for Identifying "mod.tag" Category for Each Forecast              #
###########################################################################################
id.mod.tag <- function(val.unrounded){
  modset.10 <- c(0.1,0.2,0.3,0.4,0.6,0.7,0.8,0.9)
  modset.05 <- c(0.05,0.15,0.25,0.35,0.45,0.55,0.65,0.75,0.85,0.95)
  modset.special <- c(0,0.5,1)
  modset.01 <- round(seq.int(from=0,to=100)/100,2)
  modset.01 <- modset.01[!(is.element(modset.01,modset.10))]
  modset.01 <- modset.01[!(is.element(modset.01,modset.05))]
  modset.01 <- modset.01[!(is.element(modset.01,modset.special))]

  if(val.unrounded==0){
    data.out <- "mod.0"
  }
  if(val.unrounded==0.5){
    data.out <- "mod.50"
  }
  if(val.unrounded==1){
    data.out <- "mod.1.0"
  }
  if(is.element(val.unrounded, modset.10)){
    data.out <- "mod.10"
  }
  if(is.element(val.unrounded, modset.05)){
    data.out <- "mod.05"
  }
  if(is.element(val.unrounded, modset.01)){
    data.out <- "mod.01"
  }
  
  return(data.out)
  #rm(modset.01,modset.05,modset.10,modset.special)
}


###########################################################################################
#      Function for Identifying "time.horiz.code" and "time.horiz" for Each Forecast      #
###########################################################################################
time.horiz.tag <- function(data.frame.in){
  data.frame.in <- cbind(data.frame.in, orig.order=seq.int(1,nrow(data.frame.in)))
  data.frame.in <- data.frame.in[order(data.frame.in$t.minus_days),]
  
  i <- 1
  j <- 1
  tag.vector <- vector()
  code.vector <- vector()
  
  while(i <= nrow(data.frame.in)){
    if(data.frame.in$t.minus_days[i] > 0){
      tag.vector[i] <- "post.ifp.close"
      tag.vector[i+1] <- "post.ifp.close"
      code.vector[i] <- 0
      code.vector[i+1] <- 0
    } else if(data.frame.in$t.minus_days[i] >= -14 & data.frame.in$t.minus_days[i] <= 0 & data.frame.in$val.unrounded[i] <= 0.05){
           tag.vector[i] <- "lay.up"
           tag.vector[i+1] <- "lay.up"
            code.vector[i] <- 1
            code.vector[i+1] <- 1
          }
      else if(data.frame.in$t.minus_days[i] >= -14 & data.frame.in$t.minus_days[i] <= 0 & data.frame.in$val.unrounded[i] >= 0.95){
           tag.vector[i] <- "lay.up"
           tag.vector[i+1] <- "lay.up"
           code.vector[i] <- 1
           code.vector[i+1] <- 1
          } else{
              if(j <= 259696){
                  tag.vector[i] <- "long.term"
                  tag.vector[i+1] <- "long.term"
                  code.vector[i] <- 4
                  code.vector[i+1] <- 4
                }
    
              if(j > 259696 & j <= (259696*2)){
                  tag.vector[i] <- "mid.term"
                  tag.vector[i+1] <- "mid.term"
                  code.vector[i] <- 3
                  code.vector[i+1] <- 3
                }
    
              if(j > (259696*2)){
                  tag.vector[i] <- "short.term"
                  tag.vector[i+1] <- "short.term"
                  code.vector[i] <- 2
                  code.vector[i+1] <- 2
                }
      
               j <- j+2
              }
   i <- i+2
  }
  
  
 data.out <- data.frame(time.horiz=tag.vector, time.horiz.code=code.vector)
 data.out <- data.out[order(data.frame.in$orig.order),]
 return(data.out)
  
 #rm(data.frame.in,data.out,i,j,tag.vector,code.vector)
}



###########################################################################################
#       Function for Rounding Forecasts to the Midpoint of K Equal-Sized Bins             #
###########################################################################################
Round.To <- function(data, k){
  i <- 1
  Rounded.Data <- vector()
  
  jump.size <- (1/k)
  bin.mids <- seq((jump.size/2),1,jump.size)
  
  bin.splits <- seq(0,1,jump.size)
  num.splits <- length(bin.splits) - 1
  bin.splits <- bin.splits[c(2:num.splits)]
  
  while (i <= length(data)){
    #Randomly adjust points falling at the dividing line between bins
    if(is.element(data[i], bin.splits)){
      data[i] <- data[i] + runif(1, min=-0.0000000001, max=0.0000000001)
    }
    
    diffs <- abs(bin.mids - data[i])
    Rounded.Data[i] <- bin.mids[which.min(diffs)]
    i <- i+1
  }
  
  rm(i,jump.size,bin.mids,diffs)
  return(Rounded.Data)
}


###########################################################################################
#              Function for Rounding Forecasts to the Midpoint of NIC Bins                #
###########################################################################################
Round.To_NIC <- function(data){
  i <- 1
  Rounded.Data <- vector()
  
  bin.mids <- c(0.03,0.125,0.325,0.5,0.675,0.875,0.97) ##May want to tweak this
  #bin.mids <- c(0.025,0.125,0.325,0.5,0.675,0.875,0.975)
  bin.splits <- c(0.05,0.20,0.45,0.55,0.80,0.95)
  
  while (i <= length(data)){
    #Randomly adjust points falling at the dividing line between bins
    if(is.element(data[i], bin.splits)){
      data[i] <- data[i] + runif(1, min=-0.0000000001, max=0.0000000001)
    }
    
    diffs <- abs(bin.mids - data[i])
    Rounded.Data[i] <- bin.mids[which.min(diffs)]
    i <- i+1
  }
  
  rm(i,bin.splits,bin.mids,diffs)
  return(Rounded.Data)
}


###########################################################################################
#    Function for Rounding Forecasts to the Midpoint of NIC Bins (Alternative Version)    #
###########################################################################################
Round.To_NIC.alt <- function(data){
  i <- 1
  Rounded.Data <- vector()
  
  #bin.mids <- c(0.03,0.125,0.325,0.5,0.675,0.875,0.97) ##May want to tweak this
  bin.mids <- c(0.025,0.125,0.325,0.5,0.675,0.875,0.975)
  bin.splits <- c(0.05,0.20,0.45,0.55,0.80,0.95)
  
  while (i <= length(data)){
    #Randomly adjust points falling at the dividing line between bins
    if(is.element(data[i], bin.splits)){
      data[i] <- data[i] + runif(1, min=-0.0000000001, max=0.0000000001)
    }
    
    diffs <- abs(bin.mids - data[i])
    Rounded.Data[i] <- bin.mids[which.min(diffs)]
    i <- i+1
  }
  
  rm(i,bin.splits,bin.mids,diffs)
  return(Rounded.Data)
}


###########################################################################################
#           Function for Calculating an Expected Value (Frequency Weighted Mean)          #
###########################################################################################
calc.wt.mean <- function(data.in){
  i <- 1
  weights <- vector()
  unique.vals <- unique(data.in$val.unrounded)
  
  while(i <= length(unique.vals)){
    weights[i] <- length(data.in$val.unrounded[data.in$val.unrounded==unique.vals[i]]) / nrow(data.in)
    i <- i+1
  }
  
  return(sum(unique.vals*weights))
}


###########################################################################################
#        Function for Rounding Forecasts to Bin Expected Value (Equal-Size Bins)          #
###########################################################################################
Round.To.EV <- function(data, k){
  i <- 1 #Forecast Index
  j <- 1 #Bin index
  
  bin.boundaries <- seq(from=0, to=1, 1/k)
  fcast.bin.num <- vector()
  data.out <- vector()
  
  ##Categorize data by bin##
  while(i <= length(data)){
    ##Round randomly at bin boundaries##
    if(is.element(data[i], bin.boundaries) & data[i] > 0){
      data[i] <- data[i] + runif(1, min=-0.000001, max=0.000001)
    }
    
    fcast.bin.num[i] <- min((as.integer(abs(data[i] %/% (1/k)))+1),k)
    
    ##Re-set data value as needed##
    data[i] <- round(data[i],2)
    
    i <- i+1
  }
  
  ##Re-set i##
  i <- 1
  
  ##Create data.frame##
  data <- data.frame(val.unrounded=data, bin.num=fcast.bin.num)
  
  ##Calculate Bin EVs##
  bin.evs <- as.vector(by(data, data$bin.num, calc.wt.mean))
  
  ##Round each value to bin EV##
  while(i <= nrow(data)){
    data.out[i] <- bin.evs[data$bin.num[i]]
    i <- i+1
  }
  
  data.out <- data.frame(val.ev.equi=data.out)
  
  return(data.out)
  
  #rm(data,i,j,k,bin.boundaries,fcast.bin.num,data.out,data,bin.evs,name)
}



###########################################################################################
#           Function for Rounding Forecasts to Bin Expected Value (NIC Bins)              #
###########################################################################################
Round.To.EV_NIC.7 <- function(data.frame.in){
  i <- 1 #Forecast Index
  data.out <- vector()
  
  ##Calculate Bin EVs##
  bin.evs <- as.vector(by(data.frame.in, data.frame.in$nic.bin, calc.wt.mean))
  
  ##Round each value to bin EV##
  while(i <= nrow(data.frame.in)){
    data.out[i] <- bin.evs[data.frame.in$nic.bin[i]]
    i <- i+1
  }
  
  return(data.out)
  
  #rm(data.frame.in,i,data.out,bin.evs)
}


###########################################################################################
#                        Function for Calculating Brier Scores                            #
###########################################################################################
bs <- function(dataline, value.col.index){
  if(dataline$answer_option[1]=="a"){
    pr.1 <- dataline[1,value.col.index]
    pr.0 <- 1 - dataline[1,value.col.index]
  }
  
  if(dataline$answer_option[1]=="b"){
    pr.1 <- 1 - dataline[1,value.col.index]
    pr.0 <- dataline[1,value.col.index]
  }
  
  out <- ifelse(dataline$outcome[1]=="a", (((1-pr.1)^2 + (0-pr.0)^2)/2),(((0-pr.1)^2 + (1-pr.0)^2)/2))
  return(out)
}


###########################################################################################
#                     Function for Calculating Logarithmic Scores                         #
###########################################################################################
log.score <- function(dataline, value.col.index){
  if(dataline$answer_option[1]=="a"){
    pr.1 <- dataline[1,value.col.index]
    
    ## Correct for possible infinite penalty
    if(pr.1==1){
      pr.1 <- 0.995
    }
    if(pr.1==0){
      pr.1 <- 0.005
    }
  }
  
  if(dataline$answer_option[1]=="b"){
    pr.1 <- 1 - dataline[1,value.col.index]
    ## Correct for possible infinite penalty
    if(pr.1==1){
      pr.1 <- 0.995
    }
    if(pr.1==0){
      pr.1 <- 0.005
    }
  }
  
  out <- ifelse(dataline$outcome[1]=="a", -log(pr.1),-log(1-pr.1))
  return(out)
}


###########################################################################################
#   Three-Tiered Series of Functions for Generating "Returns to Precision" Tables         #
###########################################################################################


########################### TOP LEVEL ####################################
##########################################################################

## This function represents the "top-level" of a three tiered process for generating a
## "returns to precision" table. At this level, the user can specify the scoring rule
## to be used (i.e., the "method"), the type of table to be generated (i.e., the "type"),
## and the measure of central tendency to be used in aggregation (i.e., the "stat"). Once
## specified, the "top-level" funcion calls the appropriate "middle-level" function to
## generate the specified table.
##
## ARGUMENT KEY:
##   method:
##     "bs" = brier score;
##     "log" = logarithmic scoring rule;
##     "spr" = spherical scoring rule (not used in manuscript, but included for completeness).
##   type:
##     "comp.group" = subset data according to forecaster comparison group;
##     "time.horiz" = subset data according to time horizon
##     "wep" = subset data according to equal-sized WEPs;
##     "nic" = subset data according to new NIC WEP categorization scheme."
##   stat:
##     "mean" = use the mean as an aggregator;
##     "median" = use the median as an aggregator.


gen.tab_master <- function(method=c("bs","log","spr"), type=c("comp.group", "time.horiz", "wep", "nic"), stat=c("mean","median")){  
  #Means
  if(stat=="mean"){
    if(type=="comp.group"){
      table.out <- gen.tab_comp.group(stat="mean", method=method)
    }
    if(type=="time.horiz"){
      table.out <- gen.tab_time.horiz(stat="mean", method=method)
    }
    if(type=="wep"){
      table.out <- gen.tab_wep(stat="mean", method=method)
    }
    if(type=="nic"){
      table.out <- gen.tab_nic(stat="mean", method=method)
    }
  }
  
  #Medians
  if(stat=="median"){
    if(type=="comp.group"){
      table.out <- gen.tab_comp.group(stat="median", method=method)
    }
    if(type=="time.horiz"){
      table.out <- gen.tab_time.horiz(stat="median", method=method)
    }
    if(type=="wep"){
      table.out <- gen.tab_wep(stat="median", method=method)
    }
    if(type=="nic"){
      table.out <- gen.tab_nic(stat="median", method=method)
    }
  }
  
  return(table.out)
}


############################## MID LEVEL #################################
##########################################################################

## The following four functions represent the "middle-level" of the three tiered process
## for generating a "returns to precision" table. Each of these four functions is quite
## similar to the next, but each is customized to handle the eccentricities of the particular type
## of table being generated. In each case, the general purpose of the "middle-level"
## functions is to gathe the raw materials for a "returns to precision" table (by calling the appropriate
## "low-level" function), and displaying it in an easily interpreted manner (by calling the function
## "gen.output.table").
##
## To understand these functions, a few details bear mentioning.
##
## 1.
##   In the "comp.group" function and elsewhere, comparison groups are indicated by a variable called
##   "g.tnt". This stands for "group, training/no training," and has six levels.
##      "all.fc" = all forecasters;
##      "i.t" = individuals w/training;
##      "i.nt" = individuals w/o training
##      "g.t" = non-superforecaster groups w/training;
##      "g.nt" = non-superforecaster groups w/o training;
##      "s.t" = superfoecasters, all of who were trained and worked in groups.
##
## 2.
##    In the "time.horiz" function, time periods are indicated by a variable called "time.horiz.cat",
##    which has five levels.
##      "no.lay.ups" = all forecasts that were not included in the lay-ups category (see manuscript);
##      "lay-ups" = only those forecasts that were designated lay-ups;
##      "short-term", "medium term", and "long term" = forecasts from period 1, 2, and 3, respectively
##
## 3.
##    In all four functions, the code uses shorthand to indicate observations at different levels of aggregation.
##    "iilms" = individual, ifp-level means (medians), where an ifp is a forecasting question.
##    "gilms" = group, ifp-level means (medians), where a group is determied by the type of table in question.
##    "gylms" = group, year-level means (medians). These are gilms, aggregated across ifps in a given tournamnet year.



####### Function to Generate a Comp.Group Table ##########
gen.tab_comp.group <- function(stat=c("mean", "median"), method=c("bs","log","spr")){
  i <- 1
  tags <- c("all.fc","i.nt","i.t","g.nt","g.t","s.t")
  gilm.table.list <- list()
  drop.stat.list <- list()
  gylm.table.list <- list()
  
  #Means
  if(stat=="mean"){
    while(i <= length(tags)){
      #Select Data
      sub <- fcasts
      
      if(tags[i]!="all.fc"){
        sub <- subset(sub, g.tnt==tags[i])
      }
      
      tag <- paste(method, tags[i], sep="_")
      
      #Generate Output
      temp <- gen.output(stat="mean", sub, tag, method=method)
      
      output.name <- paste("output",stat,tag, sep="_")
      assign(output.name, temp, envir=.GlobalEnv)
      
      ##Gather data for table
      gilm.table.list <- append(gilm.table.list,temp[1])
      drop.stat.list <- append(drop.stat.list,temp[2])
      gylm.table.list <- append(gylm.table.list,temp[3])
      
      ## Repeat for 2-sided test
      tag <- paste(tag, "two.sided", sep="_")
      temp <- gen.output(stat="mean", sub, tag, two.sided=TRUE, method=method)
      
      output.name <- paste("output",stat,tag, sep="_")
      assign(output.name, temp, envir=.GlobalEnv)
      
      ##Gather data for table
      gilm.table.list <- append(gilm.table.list,temp[1])
      drop.stat.list <- append(drop.stat.list,temp[2])
      gylm.table.list <- append(gylm.table.list,temp[3])
      
      
      i <- i+1
    }
  }
  
  #Medians
  if(stat=="median"){
    while(i <= length(tags)){
      #Select Data
      sub <- fcasts
      
      if(tags[i]!="all.fc"){
        sub <- subset(sub, g.tnt==tags[i])
      }
      
      tag <- paste(method, tags[i], sep="_")
      
      #Generate Output
      temp <- gen.output(stat="median", sub, tag, method=method)
      
      output.name <- paste("output",stat,tag, sep="_")
      assign(output.name, temp, envir=.GlobalEnv)
      
      ##Gather data for table
      gilm.table.list <- append(gilm.table.list,temp[1])
      drop.stat.list <- append(drop.stat.list,temp[2])
      gylm.table.list <- append(gylm.table.list,temp[3])
      
      ## Repeat for 2-sided test
      tag <- paste(tag, "two.sided", sep="_")
      temp <- gen.output(stat="median", sub, tag, two.sided=TRUE, method=method)
      
      output.name <- paste("output",stat,tag, sep="_")
      assign(output.name, temp, envir=.GlobalEnv)
      
      ##Gather data for table
      gilm.table.list <- append(gilm.table.list,temp[1])
      drop.stat.list <- append(drop.stat.list,temp[2])
      gylm.table.list <- append(gylm.table.list,temp[3])
      
      
      i <- i+1
    }
  }
  
  rm(i,tags,sub,temp,output.name)
  
  ### Display Output as a single table (function in "gen.output.table" file###
  if(stat=="mean"){
    table.out <- gen.output.table(stat="mean",gilm.table.list,drop.stat.list,gylm.table.list,bonferroni=1)
  }
  if(stat=="median"){
    table.out <- gen.output.table(stat="median",gilm.table.list,drop.stat.list,gylm.table.list,bonferroni=1)
  }
  
  return(table.out)
}



####### Function to Generate a Time.Horiz Table ##########

gen.tab_time.horiz <- function(stat=c("mean", "median"), method=c("bs","log","sphr")){
  i <- 1 #Comp.Group Index
  j <- 1 #Time.Horiz Index
  
  i.vector <- c("all.fc", "s.t")
  j.vector <- c("no.lay.ups","lay.up","short.term","mid.term","long.term")
  
  gilm.table.list <- list()
  drop.stat.list <- list()
  gylm.table.list <- list()
  
  #Means
  if(stat=="mean"){
    while(i <= length(i.vector)){
      #Select Data
      sub.i <- fcasts
      
      if(i == 2){
        sub.i <- subset(sub.i, g.tnt=="s.t")
      }
      
      while(j <= length(j.vector)){
        if(j != 1){
          sub.j <- subset(sub.i, time.horiz==j.vector[j])
        }
        if(j==1){
          sub.j <- subset(sub.i, time.horiz.code > 1)
        }
        
        tag <- paste(method,i.vector[i],j.vector[j], sep="_")
        
        #Generate Output
        temp <- gen.output(stat="mean", sub.j, tag, method=method)
        
        output.name <- paste("output",stat,tag, sep="_")
        assign(output.name, temp, envir=.GlobalEnv)
        
        ##Gather data for table
        gilm.table.list <- append(gilm.table.list,temp[1])
        drop.stat.list <- append(drop.stat.list,temp[2])
        gylm.table.list <- append(gylm.table.list,temp[3])
        
        ## Repeat for 2-sided test
        tag <- paste(tag, "two.sided", sep="_")
        temp <- gen.output(stat="mean", sub.j, tag, two.sided=TRUE, method=method)
        
        output.name <- paste("output",stat,tag, sep="_")
        assign(output.name, temp, envir=.GlobalEnv)
        
        ##Gather data for table
        gilm.table.list <- append(gilm.table.list,temp[1])
        drop.stat.list <- append(drop.stat.list,temp[2])
        gylm.table.list <- append(gylm.table.list,temp[3])
        
        j <- j+1
      }
      
      j <- 1
      i <- i+1
    }
  }
  
  #Medians
  if(stat=="median"){
    while(i <= length(i.vector)){
      #Select Data
      sub.i <- fcasts
      
      if(i == 2){
        sub.i <- subset(sub.i, g.tnt=="s.t")
      }
      
      while(j <= length(j.vector)){
        if(j != 1){
          sub.j <- subset(sub.i, time.horiz==j.vector[j])
        }
        if(j==1){
          sub.j <- subset(sub.i, time.horiz.code > 1)
        }
        
        tag <- paste(method,i.vector[i],j.vector[j], sep="_")
        
        #Generate Output
        temp <- gen.output(stat="median", sub.j, tag, method=method)
        
        output.name <- paste("output",stat,tag, sep="_")
        assign(output.name, temp, envir=.GlobalEnv)
        
        ##Gather data for table
        gilm.table.list <- append(gilm.table.list,temp[1])
        drop.stat.list <- append(drop.stat.list,temp[2])
        gylm.table.list <- append(gylm.table.list,temp[3])
        
        ## Repeat for 2-sided test
        tag <- paste(tag, "two.sided", sep="_")
        temp <- gen.output(stat="median", sub.j, tag, two.sided=TRUE, method=method)
        
        output.name <- paste("output",stat,tag, sep="_")
        assign(output.name, temp, envir=.GlobalEnv)
        
        ##Gather data for table
        gilm.table.list <- append(gilm.table.list,temp[1])
        drop.stat.list <- append(drop.stat.list,temp[2])
        gylm.table.list <- append(gylm.table.list,temp[3])
        
        
        j <- j+1
      }
      
      j <- 1
      i <- i+1
    }
  }
  
  rm(i,j,i.vector,j.vector,sub.i,sub.j,tag,sub,temp,output.name)
  
  ### Display Output as a single table (function in "gen.output.table" file###
  if(stat=="mean"){
    table.out <- gen.output.table(stat="mean",gilm.table.list,drop.stat.list,gylm.table.list,bonferroni=1)
  }
  if(stat=="median"){
    table.out <- gen.output.table(stat="median",gilm.table.list,drop.stat.list,gylm.table.list,bonferroni=1)
  }
  
  return(table.out)
}


####### Function to Generate a WEP Table ##########

gen.tab_wep <- function(stat=c("mean", "median"), method=c("bs","log","sphr")){
  i <- 1 #Comp.Group Index
  j <- 1 #WEP Index
  
  i.vector <- c("all.fc", "s.t")
  j.vector <- seq.int(from=1,to=7)
  
  gilm.table.list <- list()
  drop.stat.list <- list()
  gylm.table.list <- list()
  
  #Means
  if(stat=="mean"){
    while(i <= length(i.vector)){
      #Select Data
      sub.i <- fcasts
      
      if(i == 2){
        sub.i <- subset(sub.i, g.tnt=="s.t")
      }
      
      while(j <= length(j.vector)){
        sub.j <- subset(sub.i, wep==j.vector[j])
        
        tag <- paste(method,i.vector[i],paste("wep.", j.vector[j], sep=""), sep="_")
        
        #Generate Output
        temp <- gen.output(stat="mean", sub.j, tag, method=method)
        
        output.name <- paste("output",stat,tag, sep="_")
        assign(output.name, temp, envir=.GlobalEnv)
        
        ##Gather data for table
        gilm.table.list <- append(gilm.table.list,temp[1])
        drop.stat.list <- append(drop.stat.list,temp[2])
        gylm.table.list <- append(gylm.table.list,temp[3])
        
        ## Repeat for 2-sided test
        tag <- paste(tag, "two.sided", sep="_")
        temp <- gen.output(stat="mean", sub.j, tag, two.sided=TRUE, method=method)
        
        output.name <- paste("output",stat,tag, sep="_")
        assign(output.name, temp, envir=.GlobalEnv)
        
        ##Gather data for table
        gilm.table.list <- append(gilm.table.list,temp[1])
        drop.stat.list <- append(drop.stat.list,temp[2])
        gylm.table.list <- append(gylm.table.list,temp[3])
        
        j <- j+1
      }
      
      j <- 1
      i <- i+1
    }
  }
  
  #Medians
  if(stat=="median"){
    while(i <= length(i.vector)){
      #Select Data
      sub.i <- fcasts
      
      if(i == 2){
        sub.i <- subset(sub.i, g.tnt=="s.t")
      }
      
      while(j <= length(j.vector)){
        sub.j <- subset(sub.i, wep==j.vector[j])
        
        tag <- paste(method,i.vector[i],paste("wep.", j.vector[j], sep=""), sep="_")
        
        #Generate Output
        temp <- gen.output(stat="median", sub.j, tag, method=method)
        
        output.name <- paste("output",stat,tag, sep="_")
        assign(output.name, temp, envir=.GlobalEnv)
        
        ##Gather data for table
        gilm.table.list <- append(gilm.table.list,temp[1])
        drop.stat.list <- append(drop.stat.list,temp[2])
        gylm.table.list <- append(gylm.table.list,temp[3])
        
        ## Repeat for 2-sided test
        tag <- paste(tag, "two.sided", sep="_")
        temp <- gen.output(stat="median", sub.j, tag, two.sided=TRUE, method=method)
        
        output.name <- paste("output",stat,tag, sep="_")
        assign(output.name, temp, envir=.GlobalEnv)
        
        ##Gather data for table
        gilm.table.list <- append(gilm.table.list,temp[1])
        drop.stat.list <- append(drop.stat.list,temp[2])
        gylm.table.list <- append(gylm.table.list,temp[3])
        
        j <- j+1
      }
      
      j <- 1
      i <- i+1
    }
  }
  
  rm(i,j,i.vector,j.vector,sub.i,sub.j,tag,sub,temp,output.name)
  
  ### Display Output as a single table (function in "gen.output.table" file###
  if(stat=="mean"){
    table.out <- gen.output.table(stat="mean",gilm.table.list,drop.stat.list,gylm.table.list,bonferroni=1)
  }
  if(stat=="median"){
    table.out <- gen.output.table(stat="median",gilm.table.list,drop.stat.list,gylm.table.list,bonferroni=1)
  }
  
  return(table.out)
}



####### Function to Generate a NIC Table ##########

gen.tab_nic <- function(stat=c("mean", "median"), method=c("bs","log","sphr")){
  i <- 1 #Comp.Group Index
  j <- 1 #NIC.bin Index
  
  i.vector <- c("all.fc", "s.t")
  j.vector <- seq.int(from=1,to=7)
  
  gilm.table.list <- list()
  drop.stat.list <- list()
  gylm.table.list <- list()
  
  #Means
  if(stat=="mean"){
    while(i <= length(i.vector)){
      #Select Data
      sub.i <- fcasts
      
      if(i == 2){
        sub.i <- subset(sub.i, g.tnt=="s.t")
      }
      
      while(j <= length(j.vector)){
        sub.j <- subset(sub.i, nic.bin==j.vector[j])
        
        tag <- paste(method,i.vector[i],paste("nic.bin.", j.vector[j], sep=""), sep="_")
        
        #Generate Output
        temp <- gen.output(stat="mean", sub.j, tag, method=method)
        
        output.name <- paste("output",stat,tag, sep="_")
        assign(output.name, temp, envir=.GlobalEnv)
        
        ##Gather data for table
        gilm.table.list <- append(gilm.table.list,temp[1])
        drop.stat.list <- append(drop.stat.list,temp[2])
        gylm.table.list <- append(gylm.table.list,temp[3])
        
        ## Repeat for 2-sided test
        tag <- paste(tag, "two.sided", sep="_")
        temp <- gen.output(stat="mean", sub.j, tag, two.sided=TRUE, method=method)
        
        output.name <- paste("output",stat,tag, sep="_")
        assign(output.name, temp, envir=.GlobalEnv)
        
        ##Gather data for table
        gilm.table.list <- append(gilm.table.list,temp[1])
        drop.stat.list <- append(drop.stat.list,temp[2])
        gylm.table.list <- append(gylm.table.list,temp[3])
        
        j <- j+1
      }
      
      j <- 1
      i <- i+1
    }
  }
  
  #Medians
  if(stat=="median"){
    while(i <= length(i.vector)){
      #Select Data
      sub.i <- fcasts
      
      if(i == 2){
        sub.i <- subset(sub.i, g.tnt=="s.t")
      }
      
      while(j <= length(j.vector)){
        sub.j <- subset(sub.i, nic.bin==j.vector[j])
        
        tag <- paste(method,i.vector[i],paste("nic.bin.", j.vector[j], sep=""), sep="_")
        
        #Generate Output
        temp <- gen.output(stat="median", sub.j, tag, method=method)
        
        output.name <- paste("output",stat,tag, sep="_")
        assign(output.name, temp, envir=.GlobalEnv)
        
        ##Gather data for table
        gilm.table.list <- append(gilm.table.list,temp[1])
        drop.stat.list <- append(drop.stat.list,temp[2])
        gylm.table.list <- append(gylm.table.list,temp[3])
        
        ## Repeat for 2-sided test
        tag <- paste(tag, "two.sided", sep="_")
        temp <- gen.output(stat="median", sub.j, tag, two.sided=TRUE, method=method)
        
        output.name <- paste("output",stat,tag, sep="_")
        assign(output.name, temp, envir=.GlobalEnv)
        
        ##Gather data for table
        gilm.table.list <- append(gilm.table.list,temp[1])
        drop.stat.list <- append(drop.stat.list,temp[2])
        gylm.table.list <- append(gylm.table.list,temp[3])
        
        j <- j+1
      }
      
      j <- 1
      i <- i+1
    }
  }
  
  rm(i,j,i.vector,j.vector,sub.i,sub.j,tag,sub,temp,output.name)
  
  ### Display Output as a single table (function in "gen.output.table" file###
  if(stat=="mean"){
    table.out <- gen.output.table(stat="mean",gilm.table.list,drop.stat.list,gylm.table.list,bonferroni=1)
  }
  if(stat=="median"){
    table.out <- gen.output.table(stat="median",gilm.table.list,drop.stat.list,gylm.table.list,bonferroni=1)
  }
  
  return(table.out)
}

############################ LOW LEVEL ###################################
##########################################################################

## This function represents the "low-level" of the three-teired process for generating a
## "returns to precision table. Here, the specific statistical tests for each comparison
## are carried out, and returned to the "middle-level" for use in table construction.

gen.output <- function(stat=c("mean", "median"), source.data, tag, two.sided=FALSE, method=c("bs","log","sphr")){
  ### Generate GILM Table ###
  n <- 4
  
  if(method=="bs"){
    sub <- source.data[c(1,2,8,25:33)]
  }
  if(method=="log"){
    sub <- source.data[c(1,2,8,34:42)]
  }
  if(method=="sphr"){
    sub <- source.data[c(1,2,8,43:51)]
  }
  
  names <- colnames(sub)[4:12]
  data.skeleton <- sub[c(1,2,3)]
  gilm.table <- gilm.tab.spine(data.skeleton)
  gilm.table <- cbind(TAG=rep(tag,nrow(gilm.table)),gilm.table)
  
  while (n <= 12){
    data.in <- cbind(data.skeleton,sub[n])
    colnames(data.in)[4] <- "BS"
    new.gilms <- data.table(temp=calc.meangilms_all.fc(data.in))
    colnames(new.gilms)[1] <- names[n-3]
    
    gilm.table <- cbind(gilm.table, new.gilms)
    
    n <- n+1
  }
  
  gilm.table.name <- paste("gilm.table", tag, sep="_")
  assign(gilm.table.name, gilm.table, envir=.GlobalEnv)
  
  target.gilm.table <- gilm.table
  
  rm(n,sub,names,data.in,data.skeleton,new.gilms,gilm.table.name)
  ###
  
  
  #Generate Drop-Stat Table#
  if(method=="bs"){
    unrounded.index <- 25
  }
  if(method=="log"){
    unrounded.index <- 34
  }
  if(method=="sphr"){
    unrounded.index <- 43
  }
  
  n <- unrounded.index + 1
  sub <- source.data
  
  # Calculate gilms for raw.bs
  unrounded.in <- sub[c(1,2,8,unrounded.index)]
  colnames(unrounded.in)[4] <- "BS"
  x <- calc.meangilms_all.fc(unrounded.in)
  
  data.skeleton <- sub[c(1,2,8)]
  
  if(stat=="mean"){
    temp.dt <- data.table(UOA="temp",Group.1="raw.bs",Group.2="temp",t.stat=0,alternative="temp",p.val=0,CI.lower=0,CI.upper=0)
  }
  if(stat=="median"){
    temp.dt <- data.table(UOA="temp",Group.1="raw.bs",Group.2="temp",U.stat=0,alternative="temp",p.val=0,est.delta=0,CI.lower=0,CI.upper=0)
  }
  
  
  while (n <= (unrounded.index + 8)){
    data.in <- cbind(data.skeleton,sub[n])
    colnames(data.in)[4] <- "BS"
    y <- calc.meangilms_all.fc(data.in)
    
    if(two.sided==FALSE){
      if(stat=="mean"){
        t <- t.test(x,y,paired=TRUE,alternative="less",conf.int=TRUE)
        newrow <- data.table(UOA="Mean.GILM",Group.1="raw.bs",Group.2=colnames(sub[n]),t.stat=t$statistic,alternative=t$alternative,p.val=t$p.value,CI.lower=t$conf.int[1],CI.upper=t$conf.int[2])
      }
      
      if(stat=="median"){
        w <- wilcox.test(x,y,paired=TRUE,alternative="less",conf.int=TRUE)
        newrow <- data.table(UOA="Mean.GILM",Group.1="raw.bs",Group.2=colnames(sub[n]),U.stat=w$statistic,alternative=w$alternative,p.val=w$p.value,est.delta=w$estimate,CI.lower=w$conf.int[1],CI.upper=w$conf.int[2])
      }
    }
    
    if(two.sided==TRUE){
      if(stat=="mean"){
        t <- t.test(x,y,paired=TRUE,alternative="two.sided",conf.int=TRUE)
        newrow <- data.table(UOA="Mean.GILM",Group.1="raw.bs",Group.2=colnames(sub[n]),t.stat=t$statistic,alternative=t$alternative,p.val=t$p.value,CI.lower=t$conf.int[1],CI.upper=t$conf.int[2])
      }
      
      if(stat=="median"){
        w <- wilcox.test(x,y,paired=TRUE,alternative="two.sided",conf.int=TRUE)
        newrow <- data.table(UOA="Mean.GILM",Group.1="raw.bs",Group.2=colnames(sub[n]),U.stat=w$statistic,alternative=w$alternative,p.val=w$p.value,est.delta=w$estimate,CI.lower=w$conf.int[1],CI.upper=w$conf.int[2])
      }
    }
    
    temp.dt <- rbind(temp.dt,newrow)
    n <- n+1
  }
  
  temp.dt <- temp.dt[-1]
  
  if(stat=="mean"){
    drop.table.name <- paste("drop.t", tag, sep="_")
  }
  if(stat=="median"){
    drop.table.name <- paste("drop.w", tag, sep="_")
  }
  
  assign(drop.table.name, temp.dt, envir=.GlobalEnv)
  
  target.drop.table <- temp.dt
  
  rm(n,x,y,w,t,data.in,unrounded.in,sub,newrow,data.skeleton,temp.dt,drop.table.name)
  ###
  
  
  #Calculate GYLMS
  n <- 4
  sub <- source.data[c(1,2,8,unrounded.index:(unrounded.index+8))]
  
  names <- colnames(sub)[4:12]
  data.skeleton <- sub[c(1,2,3)]
  gylms <- vector()
  
  while (n <= 12){
    data.in <- cbind(data.skeleton,sub[n])
    colnames(data.in)[4] <- "BS"
    gylms[n-3] <- calc.hybridgylms_all.fc(data.in)
    
    n <- n+1
  }
  
  target.gylm.table <- data.table(subset=names, GYLM=gylms)
  gylm.table.name <- paste("hybridgylms", tag, sep="_")
  assign(gylm.table.name, target.gylm.table, envir=.GlobalEnv)
  
  rm(gylm.table.name,data.in,n,sub,names,data.skeleton,gylms)
  ###
  
  data.out <- list(target.gilm.table,target.drop.table,target.gylm.table)
  rm(target.gilm.table,target.drop.table,target.gylm.table)
  #rm(stat,source.data,tag,WEP)
  
  return(data.out)
}



################## CREATE FINAL TABLE ####################################
##########################################################################

## This function recombines and reorganizes the statistical results from the "low-level"
## into an easily interpreted output table. This table is constructed in the "middle-level"
## function, and passed back-up to the "high-level" as the final output of the gen.tab_master function.

gen.output.table <- function(stat, gilms.in, drop.stat.in, gylms.in, bonferroni = 1){
  table.out <- data.frame(stat="temp",subset="temp",mean.raw.score=0,median.raw.score=0,k.7="temp",k.5="temp",k.3="temp",k.2="temp",nic.mid="temp",nic.alt.mid="temp",ev.equi="temp",ev.nic="temp")
  
  i <- 1 #Comp.Group Index
  
  if(stat=="mean"){
    while(i <= length(gilms.in)){
      target.gilms <- as.data.frame(gilms.in[i])
      target.drop.stat <- as.data.frame(drop.stat.in[i])
      target.gylms <- as.data.frame(gylms.in[i])
      
      val.7 <- round(mean((target.gilms[,3] - target.gilms[,7]) / target.gilms[,3])*100,2)
      sig.7 <- target.drop.stat$p.val[4]*bonferroni
      stars.7 <- add.stars(sig.7)
      cell.7 <- paste(val.7, "%", sep="")
      cell.7 <- paste(cell.7, stars.7, sep="")
      
      val.5 <- round(mean((target.gilms[,3] - target.gilms[,6]) / target.gilms[,3])*100,2)
      sig.5 <- target.drop.stat$p.val[3]*bonferroni
      stars.5 <- add.stars(sig.5)
      cell.5 <- paste(val.5, "%", sep="")
      cell.5 <- paste(cell.5, stars.5, sep="")
      
      val.3 <- round(mean((target.gilms[,3] - target.gilms[,5]) / target.gilms[,3])*100,2)
      sig.3 <- target.drop.stat$p.val[2]*bonferroni
      stars.3 <- add.stars(sig.3)
      cell.3 <- paste(val.3, "%", sep="")
      cell.3 <- paste(cell.3, stars.3, sep="")
      
      val.2 <- round(mean((target.gilms[,3] - target.gilms[,4]) / target.gilms[,3])*100,2)
      sig.2 <- target.drop.stat$p.val[1]*bonferroni
      stars.2 <- add.stars(sig.2)
      cell.2 <- paste(val.2, "%", sep="")
      cell.2 <- paste(cell.2, stars.2, sep="")
      
      val.nm <- round(mean((target.gilms[,3] - target.gilms[,8]) / target.gilms[,3])*100,2)
      sig.nm <- target.drop.stat$p.val[5]*bonferroni
      stars.nm <- add.stars(sig.nm)
      cell.nm <- paste(val.nm, "%", sep="")
      cell.nm <- paste(cell.nm, stars.nm, sep="")
      
      val.nam <- round(mean((target.gilms[,3] - target.gilms[,9]) / target.gilms[,3])*100,2)
      sig.nam <- target.drop.stat$p.val[6]*bonferroni
      stars.nam <- add.stars(sig.nam)
      cell.nam <- paste(val.nam, "%", sep="")
      cell.nam <- paste(cell.nam, stars.nam, sep="")
      
      val.eve <- round(mean((target.gilms[,3] - target.gilms[,10]) / target.gilms[,3])*100,2)
      sig.eve <- target.drop.stat$p.val[7]*bonferroni
      stars.eve <- add.stars(sig.eve)
      cell.eve <- paste(val.eve, "%", sep="")
      cell.eve <- paste(cell.eve, stars.eve, sep="")
      
      val.evn <- round(mean((target.gilms[,3] - target.gilms[,11]) / target.gilms[,3])*100,2)
      sig.evn <- target.drop.stat$p.val[8]*bonferroni
      stars.evn <- add.stars(sig.evn)
      cell.evn <- paste(val.evn, "%", sep="")
      cell.evn <- paste(cell.evn, stars.evn, sep="")
      
      
      newrow <- data.frame(stat=stat,subset=target.gilms$TAG[1],mean.raw.score=mean(target.gilms[,3]),median.raw.score=median(target.gilms[,3]),k.7=cell.7,k.5=cell.5,k.3=cell.3,k.2=cell.2,nic.mid=cell.nm,nic.alt.mid=cell.nam,ev.equi=cell.eve,ev.nic=cell.evn)
      table.out <- rbind(table.out,newrow)  
      
      i <- i+1
    }
  }
  
  if(stat=="median"){
    while(i <= length(gilms.in)){
      target.gilms <- as.data.frame(gilms.in[i])
      target.drop.stat <- as.data.frame(drop.stat.in[i])
      target.gylms <- as.data.frame(gylms.in[i])
      
      val.7 <- round(median((target.gilms[,3] - target.gilms[,7]) / target.gilms[,3])*100,2)
      sig.7 <- target.drop.stat$p.val[4]*bonferroni
      stars.7 <- add.stars(sig.7)
      cell.7 <- paste(val.7, "%", sep="")
      cell.7 <- paste(cell.7, stars.7, sep="")
      
      val.5 <- round(median((target.gilms[,3] - target.gilms[,6]) / target.gilms[,3])*100,2)
      sig.5 <- target.drop.stat$p.val[3]*bonferroni
      stars.5 <- add.stars(sig.5)
      cell.5 <- paste(val.5, "%", sep="")
      cell.5 <- paste(cell.5, stars.5, sep="")
      
      val.3 <- round(median((target.gilms[,3] - target.gilms[,5]) / target.gilms[,3])*100,2)
      sig.3 <- target.drop.stat$p.val[2]*bonferroni
      stars.3 <- add.stars(sig.3)
      cell.3 <- paste(val.3, "%", sep="")
      cell.3 <- paste(cell.3, stars.3, sep="")
      
      val.2 <- round(median((target.gilms[,3] - target.gilms[,4]) / target.gilms[,3])*100,2)
      sig.2 <- target.drop.stat$p.val[1]*bonferroni
      stars.2 <- add.stars(sig.2)
      cell.2 <- paste(val.2, "%", sep="")
      cell.2 <- paste(cell.2, stars.2, sep="")
      
      val.nm <- round(median((target.gilms[,3] - target.gilms[,8]) / target.gilms[,3])*100,2)
      sig.nm <- target.drop.stat$p.val[5]*bonferroni
      stars.nm <- add.stars(sig.nm)
      cell.nm <- paste(val.nm, "%", sep="")
      cell.nm <- paste(cell.nm, stars.nm, sep="")
      
      val.nam <- round(median((target.gilms[,3] - target.gilms[,9]) / target.gilms[,3])*100,2)
      sig.nam <- target.drop.stat$p.val[6]*bonferroni
      stars.nam <- add.stars(sig.nam)
      cell.nam <- paste(val.nam, "%", sep="")
      cell.nam <- paste(cell.nam, stars.nam, sep="")
      
      val.eve <- round(median((target.gilms[,3] - target.gilms[,10]) / target.gilms[,3])*100,2)
      sig.eve <- target.drop.stat$p.val[7]*bonferroni
      stars.eve <- add.stars(sig.eve)
      cell.eve <- paste(val.eve, "%", sep="")
      cell.eve <- paste(cell.eve, stars.eve, sep="")
      
      val.evn <- round(median((target.gilms[,3] - target.gilms[,11]) / target.gilms[,3])*100,2)
      sig.evn <- target.drop.stat$p.val[8]*bonferroni
      stars.evn <- add.stars(sig.evn)
      cell.evn <- paste(val.evn, "%", sep="")
      cell.evn <- paste(cell.evn, stars.evn, sep="")
      
      
      newrow <- data.frame(stat=stat,subset=target.gilms$TAG[1],mean.raw.score=mean(target.gilms[,3]),median.raw.score=median(target.gilms[,3]),k.7=cell.7,k.5=cell.5,k.3=cell.3,k.2=cell.2,nic.mid=cell.nm,nic.alt.mid=cell.nam,ev.equi=cell.eve,ev.nic=cell.evn)
      table.out <- rbind(table.out,newrow)  
      
      i <- i+1
    }
  }
  
  table.out <- table.out[-1,]
  
  return(table.out)
  #rm(stat,gilms.in,drop.stat.in,gylms.in,bonferroni,table.out,i,val.7,sig.7,stars.7,cell.7,val.5,sig.5,stars.5,cell.5,val.3,sig.3,stars.3,cell.3,val.2,sig.2,stars.2,cell.2,val.nm,sig.nm,stars.nm,cell.nm,val.nam,sig.nam,stars.nam,cell.nam,val.eve,sig.eve,stars.eve,cell.eve,val.evn,sig.evn,stars.evn,cell.evn,newrow)
}


###########################################################################################
#                                Supporting Functions                                     #
###########################################################################################

## The following functions are used internally to several of the other functions in this document


### gilm.tab.spine ####
gilm.tab.spine <- function(data.in){
  require(data.table)
  ifps.list <- as.factor(unique(data.in$ifp_id)) #Active IFPs in tournament year
  ifps.list <- ifps.list[order(ifps.list)] #Sort IFPS
  
  return(data.table(ifp.id=ifps.list))
}
####


#### add.stars ####
add.stars <- function(p.val){
  if(p.val > 0.05){stars <- c("")}
  if(p.val <= 0.05){stars <- c(" *")}
  if(p.val < 0.01){stars <- c(" **")}
  if(p.val < 0.001){stars <- c(" ***")}
  
  return(noquote(stars))
}

###### Calculate HYBRIDGYLMs for ALL FORECASTERS
# requires: ct, user_id, ifp_id, value
calc.hybridgylms_all.fc <- function(data.in){
  i <- 1 #IFP index
  
  gilms <- vector()
  
  ifps.list <- as.factor(unique(data.in$ifp_id)) #Active IFPs in tournament year
  ifps.list <- ifps.list[order(ifps.list)] #Sort IFPS
  
  while (i <= length(ifps.list)){
    sub.ifp <- subset(data.in, ifp_id==ifps.list[i])
    sub.ifp$user_id <- as.factor(as.numeric(as.character(sub.ifp$user_id))) #re-set number of factor levels
    sub.ifp <- sub.ifp[order(sub.ifp$user_id),]
    
    #Active users on IFP
    user.list <- as.factor(unique(as.numeric(as.character(sub.ifp$user_id))))
    user.list <- user.list[order(user.list)]
    
    ifp.iilms <- data.table(user_id=user.list)
    ifp.iilms <- cbind(ifp.iilms, ifp.mean=as.vector(tapply(sub.ifp$BS, sub.ifp$user_id, mean)))
    
    ## GILMS
    gilms[i] <- mean(ifp.iilms$ifp.mean) 
    
    i <- i+1
  }
  
  gylm <- median(gilms)
  
  return(gylm)
}

#For removal of vars after walk-through testing
#rm(i,gilms,ifps.list,sub.ifp,user.list,ifp.iilms,gylm)


calc.meangilms_all.fc <- function(data.in){
  i <- 1 #IFP index
  
  gilms <- vector()
  
  ifps.list <- as.factor(unique(data.in$ifp_id)) #Active IFPs in tournament year
  ifps.list <- ifps.list[order(ifps.list)] #Sort IFPS
  
  while (i <= length(ifps.list)){
    sub.ifp <- subset(data.in, ifp_id==ifps.list[i])
    sub.ifp$user_id <- as.factor(as.numeric(as.character(sub.ifp$user_id))) #re-set number of factor levels
    sub.ifp <- sub.ifp[order(sub.ifp$user_id),]
    
    #Active users on IFP
    user.list <- as.factor(unique(as.numeric(as.character(sub.ifp$user_id))))
    user.list <- user.list[order(user.list)]
    
    ifp.iilms <- data.table(user_id=user.list)
    ifp.iilms <- cbind(ifp.iilms, ifp.mean=as.vector(tapply(sub.ifp$BS, sub.ifp$user_id, mean)))
    
    ## GILMS
    gilms[i] <- mean(ifp.iilms$ifp.mean) 
    
    i <- i+1
  }
  
  return(gilms)
}

#For removal of vars after walk-through testing
#rm(i,gilms,ifps.list,sub.ifp,user.list,ifp.iilms)
###